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Abstract. While abstract interpretation is not theoretically restricted 
to specific kinds of properties, it is, in practice, mainly developed to 
compute linear over-approximations of reachable sets, aka. the collecting 
semantics of the program. The verification of user-provided properties is 
not easily compatible with the usual forward fixpoint computation using 
numerical abstract domains. 

We propose here to rely on sums-of-squares programming to characterize 
a property-driven polynomial invariant. This invariant generation can be 
guided by either boundedness, or in contrary, a given zone of the state 
space to avoid. 

While the target property is not necessarily inductive with respect to 
the program semantics, our method identifies a stronger inductive poly¬ 
nomial invariant using numerical optimization. Our method applies to a 
wide set of programs: a main while loop composed of a disjunction (if- 
then-else) of polynomial updates e.g. piecewise polynomial controllers. 
It has been evaluated on various programs. 


1 Introduction 

With the increased need for confidence in software, it becomes more than ever im¬ 
portant to provide means to support the verification of specification of software. 
Among the various formal verification methods to support these analysis, a first 
line of approaches, such as deductive methods or SMT-based model checking, 
provide rich languages to support the expression of the specification and then try 
to discharge the associate proof obligation using automatic solvers. The current 
state of the art of these solvers is able to manipulate satisfiability problems over 
linear arithmetics or restricted fragments of non linear arithmetics. Another line 
of approaches, such as static analysis also known as abstract interpretation, re¬ 
stricts, a priori, the kind of properties considered during the computation: these 
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methods typically perform interval arithmetic analysis or rely on convex poly- 
hedra computations. In practice this second line of work seems more capable of 
manipulating and generating numerical invariants through the computation of 
inductive invariants, while the first line of approaches hardly synthesize these 
required invariants through satisfiability checks. 

However, when it comes to more than linear properties, the state of the art 
is not well developed. In the early 2000s, ellipsoid analyses [Fer041 . similar to 
restricted cases of Lyapunov functions, were designed to support the study of a 
family of Airbus controllers. This exciting result was used to provide the analysis 
of absence of runtime errors but could hardly be adapted to handle more general 
user provided specifications for polynomial programs. 

However proving polynomial inequalities is NP-hard and boils down to show 
that the infimum of a given polynomial is nonnegative. Still, one can obtain lower 
bounds of such infima by decomposing certain nonnegative polynomials into 
sums-of-squares (SOS). This actually leads to solve hierarchies of semidefinite 
relaxations, introduced by Lasserre in [Las M- Recent advances in semidefinite 
programming allowed to extensively apply these relaxations to various fields, 
including parametric polynomial optimization, optimal control, combinatorial 
optimization, etc. (see e.g. jPar03ILau09| for more details). 

While these approaches were mentioned a decade ago in |Cou05| and mainly 
applied to termination analysis, they hardly made their way through the software 
verification community to address more general properties. 


Contributions. Our contribution allows to analyze high level properties defined 
as a sublevel set of polynomials functions, i.e. basic semialgebraic sets. This 
class of properties is rather large: it ranges from boundedness properties to the 
definition of a bad region of the state space to avoid. While these properties, 
when they hold, are meant to be invariant, i.e. they hold in each reachable state, 
they are not necessarily inductive. Our approach rely on the computation of 
a stronger inductive property using SOS programming. This stronger property 
is proved inductive on the complete system and, by construction, implies the 
target property specified by the user. We develop our analysis on discrete-time 
piecewise polynomial systems, capturing a wide class of critical programs, as 
typically found in current embedded systems such as aircrafts. 


Organization of the paper. The paper is organized as follows. In Section [2j we 
present the programs that we want to analyze and their representation as piece- 
wise polynomial discrete-time systems. Next, we recall in Section [3] the collect¬ 
ing semantics that we use and introduce the polynomial optimization problem 
providing inductive invariants based on target polynomial properties. Section [4] 
contains the main contribution of the paper, namely how to compute effectively 
such invariants with SOS programming. Practical computation examples are 
provided in Section [5] Finally, we explain in Section [6] how to derive template 
bases from generated invariants. 










2 Polynomial programs and piecewise polynomial 
discrete-time systems 

In this section, we describe the programs which are considered in this paper 
and we explain how to analyze them through their representation as piecewise 
polynomial discrete-time dynamical systems. 

We focus on programs composed of a single loop with a possibly complicated 
switch-case type loop body. Moreover we suppose without loss of generality that 
the analyzed programs are written in Static Single Assignment (SSA) form, that 
is each variable is initialized at most once. 

Definitions. We recall that a function / from R d to R is a polynomial if and only 
if there exists k £ N, a family {c a \ a = (aq,..., ad) £ N d , |a| = a\ + ... + ad < 
k} such that for all x £ R d , /( x) = X]| Q i</c c a.x°f 1 ... x°f d . By extension a function 
/ : R d K > R d is a polynomial if and only if all its coordinate functions are 
polynomials. Let R[x] stands for the set of d-variate polynomials. 

In this paper, we consider assignments of variables using only parallel polyno¬ 
mial assignments {x\, ..., Xd) = T(x i,..., x d) where (aq, ..., Xd) is the vector of 
the program variables. Tests are either weak polynomial inequalities r(x i,..., Xd) < 
0 or strict polynomial inequalities r(x i,..., Xd) < 0. We assume that assignments 
are polynomials from R d to R d and test functions are polynomials from R d to 
R. In the program syntax, the notation <C will be either <= or <. The form of 
the analyzed program is described in Figure [T] 
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Fig. 1 . One-loop programs with nested conditional branches 


A set C C R d is said to be a basic semialgebraic set if there exist gi,, g m £ 
R[x] such that C = {x £ R d | gj{x) <C 0, V j = 1,..., m}, where <C is used to 
encode either a strict or a weak inequality. 

As depicted in Figure [l] an update T* : R d —> R d of the *-th condition branch 
is executed if and only if the conjunction of tests r® (x) <C 0 holds. In other words, 
the variable x is updated by T l (x) if the current value of x belongs to the basic 
semialgebraic set 


X 7 := {x £ R d | Vj = 1,..., rii, r l j(x) <C 0} . 


(1) 



Piecewise Polynomial Systems. Consequently, we interpret programs as con¬ 
strained piecewise polynomial discrete-time dynamical systems (PPS for short). 
The term piecewise means that there exists a partition {X % ,i £ X } of such 
that for all i £ I, the dynamics of the system is represented by the following 
relation, for k £ N: 

if x k £ X 1 n X°, x k+1 = T\x k ). (2) 

We assume that X is finite and that the initial condition Xq belongs to some 
compact basic semialgebraic set X ln . For the program, X m is the set where 
the variables are supposed to be initialized in. Since the test entry for the loop 
condition can be nontrivial, we add the term constrained and X° denotes the 
set representing the conjunctions of tests for the loop condition. The iterates of 
the PPS are constrained to live in X°: if for some step k £ N, x k £ X° then the 
PPS is stopped at this iterate with the terminal value x k - 

We define a partition as a family of nonempty sets such that: 

U W = R d , Vi, j £ X , i ± j, X 1 n X* ± 0. (3) 

iex 

From Equation ([3]), for all k £ N* there exists a unique i £ X such that x k £ X 1 . 
A set X 1 can contain both strict and weak polynomial inequalities and character¬ 
izes the set of the m conjunctions of tests polynomials r). Let r l = (r\,... ,r z n .) 
stands for the vector of tests functions associated to the set X 1 . We suppose that 
the basic semialgebraic sets X ln and X° also admits the representation given by 
Equation 0 and we denote by r° the vector of tests polynomials (rj,... , r° o ) 
and by r ln the vector of test polynomials (r) n ,... ,r£) ). To sum up, we give a 
formal definition of PPS. 

Definition 1 (PPS). A constrained polynomial piecewise discrete-time dynam¬ 
ical system (PPS) is the quadruple (X ln ,X°, A,£) with: 

— X m C is the compact basic semialgebraic set of the possible initial con¬ 
ditions; 

— X° C R d is the basic semialgebraic set where the state variable lives; 

— X := {X l ,i £ X} is a partition as defined in Equation @1 , 

— £ := {T l ,i £ X} is the family of the polynomials from to w.r.t. the 
partition X satisfying Equation 0. 

From now on, we associate a PPS representation to each program of the form 
described at Figure [T] Since a program admits several PPS representations, we 
choose one of them, but this arbitrary choice does not change the results provided 
in this paper. In the sequel, we will often refer to the running example described 
in Example [I] 

Example 1 (Running example). The program below involves four variables and 
contains an infinite loop with a conditional branch in the loop body. The update 
of each branch is polynomial. The parameters Cij (resp. dij) are given parameters. 
During the analysis, we only keep the variables X\ and X 2 since oldx\ and oldx 2 
are just memories. 


xi,x 2 G [a\ ,a 2 ] x [6 i,6 2 ] ; 
oldx 1 = X \; 
oldx 2 = £2; 
while (-1 <= 0 ){ 
oldx 1 = £1; 
oMa;2 = X 2 ; 

case : oldx i ~2 + oldx 2^2 <= 1 : 

£1 = cn * oldx \^2 + cn * oldx 2 3 ; 

X2 = C21 * oldx ^3 + C22 * oldx 2~2; 

case : —oldxi ''2 — oldx ^2 < —1 

aq = dii * oldx i ^3 + di2 * oldx^ 2 \ 

X2 = <^21 * Oldx 1~2 + C?22 * oldx2^2\ 

} 

} 

The associated PPS corresponds to the quadruple (X m , X 0 , {X 1 , X 2 }, {T 1 , T 2 }), 
where the set of initial conditions is: 


X m — [ai,a 2 ] x [61,62], 

the system is not globally constrained, i.e. the set X° in which the variable 
x = (aq, X 2 ) lies is: 


x u = : 


the partition verifying Equation © is: 


X 1 = {a: G 


< 1 }, X 2 = {x G R 2 | —x\ — x 2 < — 1 } , 


and the polynomials relative to the partition {X 1 ,^ 2 } are: 

T\x) = ( Cll i + Cl2 i) and 7» = (+ 

\C2ixf + C22X2 J \d 2 ix( + d 2 2xi J 


3 Program invariants as sublevel sets 

The main goal of the paper is to decide automatically if a given property holds 
for the analyzed program, i.e. for all its reachable states. We are interested in 
numerical properties and more precisely in properties on the values taken by the 
d-uplet of the variables of the program. Hence, in our point-of-view, a property 
is just the membership of some set P C K d . In particular, we study properties 
which are valid after an arbitrary number of loop iterates. Such properties are 
called loop invariants of the program. Formally, we use the PPS representation 
of a given program and we say that P is a loop invariant of this program if: 

V k G N, x k G P , 

where Xk is defined at Equation © as the state variable at step k G N of the PPS 
representation of the program. Our approach addresses any property expressible 
as a polynomial level set property. This section defines formally these notions 
and develop our approach: synthesize a property-driven inductive invariant. 


3.1 Collecting Semantics as postfixpoint characterization 

Now, let us consider a program of the form described in Figure |T] and let us 
denote by S the PPS representation of this program. The set 91 of reachable 
values is the set of all possible values taken by the state variable along the 
running of S. We define 91 as follows: 

*=lK x o(* in ) ( 4 ) 

fee n 

where Tj x0 is the restriction of T on X° and Xj y0 is not defined outside X°. To 
prove that a set P is a loop invariant of the program is equivalent to prove that 
DiC. P. We can rewrite 91 inductively: 

91 = x in u (J t (91 n X 1 n x°) . ( 5 ) 

iex 

Let us denote by p(R d ) the set of subsets of R d and introduce the map F : 
p(R d ) —► p(R d ) defined by: 

F(C) = x in u |J T (c n x* n x°) (6) 

iei 

We equip p(R d ) with the partial order of inclusion. The infimum is understood 
in this sense i.e. as the greatest lower bound with respect to this order. The 
smallest fixed point problem is: 

inf {Cg p(R d ) | C = F(C)} . 

It is well-known from Tarski’s theorem that the solution of this problem exists, 
is unique and in this case, it corresponds to 91. Tarski’s theorem also states that 
91 is the smallest solution of the following Problem: 

inf {C G p(R d ) | F{C) C C} . 

Note also that the map F corresponds to a standard transfer function (or 
collecting semantics functional) applied to the PPS representation of a program. 
We refer the reader to [CC77j for a seminal presentation of this approach. 

To prove that a subset P is a loop invariant, it suffices to show that P satisfies 
F(P) C P. In this case, such P is called inductive invariant. 

3.2 Considered properties: sublevel properties P K , a 

In this paper, we consider special properties: those that are encoded with sublevel 
sets of a given polynomial function. 

Definition 2 (Sublevel property). Given a polynomial function k £ R[x] and 
aelll {Too}, we define the sublevel property V K>a as follows: 

Pk.u ■= {x G R d | k(x) <C a} . 

where <C denotes < when a G 1 and denotes < for + 00 . The expression k(x) < 
+00 expresses the boundedness of k(x) without providing a specific bound a. 




Example 2 (Sublevel property examples). 

Boundedness. When one wants to bound the reachable values of a system, we 
can try to bound the Z 2 -norm of the system: ’Py.|| 2 i00 with k(x) = ||cc|||. The use 
of a = oo does not impose any bound on k(x). 

Safe set. Similarly, it is possible to check whether a specific bound is matched. 
Either globally using the Z 2 -norm and a specific a: P||.|| 2 )Q: , or bounding the 
reachable values of each variable: P KijQi with k* : x K > Xi and at € R. 

Avoiding bad regions. If the bad region can be encoded as a sublevel property 
k(x) < 0 then its negation —k(x) < 0 characterize the avoidance of that bad 
zone. Eg. if one wants to prove that the square norm of the program variables 
is always greater than 1, then we can consider the property V K ^ a with k(x) = 
1 — ||a;||| and a = 0. 

A sublevel property is called sublevel invariant when this property is a 
loop invariant. This turns out to be difficult to prove loop invariant properties 
while considering directly 91, thus we propose to find a more tractable over¬ 
approximation of 91 for which such properties hold. 


3.3 Approach: compute a "P^a-driven inductive invariant P 

In this subsection, we explain how to compute a d -variate polynomial p e l[i] 
and a bound w £ M, such that the polynomial sublevel sets P := {x £ M. d \ 
p(x) < 0} and V K , W satisfy: 


n?C PCP c V 

1 _ 1 2= t K,w — ' K,Ot 


( 7 ) 


The first (from the left) inclusion forces P to be valid for the whole reachable 
values set. The second inclusion constraints all elements of P to satisfy the 
given sublevel property for a certain bound w. The last inclusion requires that 
the bound w is smaller than the desired level a. When a = oo, any bound w 
ensures the sublevel property. 

Now, we derive sufficient conditions on p and w to satisfy Equation ([7]). We 
decompose the problem in two parts. To satisfy the first inclusion, i.e. ensure 
that P is a loop invariant, it suffices to guarantee that F(P) C P, namely that 
P is an inductive invariant. Using Equation ([5]), P is an inductive invariant if 
and only if: 

x in u (J T (p n x i n x°) c p, 

iei 


or equivalently: 


f X in CP, 

\v* ei, r(pn x l ni°)cp. 


( 8 ) 


Thus, we obtain: 


p(x) < 0 , Wx£ X in , 

Vi e P, p{r i (x)) < o, ViePnrni 0 . 


(9) 


Now, we are interested in the second and third inclusions at Equation ([7| that 
is the sublevel property satisfaction. The condition P C V K _ W C V Ka can be 
formulated as follows: 


k(x) < w < a , \/x £ P. 


( 10 ) 


We recall that we have supposed that P is written as {x £ | p(x) < 0} where 

p £ R[x]. Finally, we provide sufficient conditions to satisfy both |9]) and ( [TO} . 
Consider the following optimization problem: 


! inf P GR[x],^eR 5 

s.t. p(x) < 0 , 

V i £l,p(T i (x)) < p{x ), 
k(x) < w + p(x), 


Vx £ X in , 

Vx e n x °, 
Vx <E R d . 


( 11 ) 


We remark that a is not present in Problem (111. Indeed, since we mini¬ 
mize w, either there exists a feasible w such that w < a and we can exploit 
this solution or such w is not available and we cannot conclude. However, from 
Problem (111, we can extract (p, w) and in the case where the optimal bound w 
is greater than a , we could use this solution with another method such as policy 
iteration [AGM15] . 


Lemma 1. Let ( p,w ) be any feasible solution of Problem ( |11[) with w < a or 
w < 00 in the case of a = 00 . Then ( p,w ) satisfies both (|9| and (10) with 
P := {x £R. d | p(x) < 0}. Finally, P and V KlW satisfy Equation (13- 


In practice, we rely on sum-of-squares programming to solve a relaxed version 


of Problem (111. 


4 Sums-of-Squares Programming for Invariant 
Generation 

We first recall some basic background about sums-of-squares certificates for poly¬ 
nomial optimization. Let K[a;] 2 m stands for the set of polynomials of degree at 
most 2m and E[x\ C R[x] be the cone of sums-of-squares (SOS) polynomials, 
that is E\x\ := { qf , with q t £ R[ar] }. Our work will use the simple fact that 
for all p £ E[x\, then p(x) > 0 for all x £ i.e. E[x\ is a restriction of the 
set of the nonnegative polynomials. For q £ R[x] 2 m > finding a SOS decomposi¬ 
tion q = qf valid over is equivalent to solve the following matrix linear 
feasibility problem: 


q{x) = bm{x) T Q b m {x) , Vx £ R d , (12) 

where b m (x) := (1, x ±,..., Xd, xf, X 1 X 2 , ■ ■ ■, x™) (the vector of all monomials 
in x up to degree m) and Q being a semidefinite positive matrix (i.e. all the 
eigenvalues of Q are nonnegative). The size of Q (as well as the length of b m ) is 

CD- 









Example 3. consider the bi-variate polynomial q(x) := l + x\ — 2x\X2 + x\. With 
b\(x) = (\,Xi,X 2 ), one looks for a semidefinite positive matrix Q such that the 
polynomial equality q(x) = bi(x) T Qb 1 (x) holds for all x G R 2 . The matrix 



satisfies this equality and has three nonnegative eigenvalues, which are 0,1, and 
2, respectively associated to the three eigenvectors eo := (0,1,1) T , e\ := (1,0,0) T 
and e 2 := (0,1,— 1) T . Defining the matrices L := (ei e 2 eo) = and 

D = ^ 02 oj, one obtains the decomposition Q = L T D L and the equality 

q(x) = {Lbi{x)) T D (Lbi{x)) = a(x) = 1 + (x\ — X 2) 2 , for all x G R 2 . The 
polynomial a is called a SOS certificate and guarantees that q is nonnegative. 


In practice, one can solve the general problem (121 by using semidefinite pro¬ 


gramming (SDP) solvers (e.g. Mosek i AAOO !. SDPA YFN + 10 ). For more de¬ 
tails about SDP, we refer the interested reader to [VB941 . 


One way to strengthen the three nonnegativity constraints of Problem (111 is 


to consider the following hierarchy of SOS programs, parametrized by the integer 
m representing the half of the degree of p: 


inf 

K [a? ] 2 m, "tu £ K 

S.t. 


- p = (To - ^2 CTl ' r 


O'3 


3=1 


Vi Gl, p~poT l = cr l ~Y p)r) - Y l)r 


0 

3 > 


3=1 


3=1 


W + p — K = if , 


V j = 1,... ,n in , crj G E[x\ , deg(<Tjr‘ n ) < 2m , 

(To G E[x] , deg(cr 0 ) < 2 m , 

V i Gl , a 1 e E[x\ , deg(cP) < 2m deg T* , 

Vi el , V j = 1,... ,Tii , pit e S[x] , deg(/i*-r*-) < 2m deg T* , 
Vi el , V j = 1,... ,n 0 , Ye E[x] , deg( 7 *r°) < 2mdegT I , 
if e E[x\ , deg(0) < 2m . 


(13) 


Proposition 1. For a given m G N, let (p m ,w m ) be any feasible solution of 
Problem ©. Then (pm,w m ) is also a feasible solution of Problem ( |11[ ). More¬ 
over, if w m < a then both P m := {1 £ R <i Pm(x) < 0} and V K>Wm satisfy 
Equation ([7|. 











Proof. The feasible solution (p m ,Wm) is associated with SOS certificates ensur¬ 


ing that the three equality constraints of Problem (131 hold: {cr 0l cr,} is associated 
to the first one, {cr*,/i*-, 7 *} is associated to the second one and if is associated 
to the third one. The first equality constraint, namely 


- Pm{x ) = a 0 (x) - 

j=1 


>( x ) r j 1 ( x ), 


Vx £ 


implies that Vx £ X ln ,p. m (x) < 0. Similarly, one has Vi £ I, Vx £ X 1 fl 
X°,p m (T l (x)) < p m {x) and Vx £ K d ,/c(x) < + p m (x). Then {p mi w m ) 


is a feasible solution of Problem (11). The second statement comes directly from 
Lemma [T] 


Computational considerations. Define t := maxjdegT*,?’ £ X\. At step m of 
this hierarchy, the number of SDP variables is proportional to ( + j mt ) and the 
number of SDP constraints is proportional to ( d+ J nt ) ■ Thus, one expects tractable 
approximations when the number d of variables (resp. the degree 2 m of the 


template p) is small. However, one can handle bigger instances of Problem (13) by 


taking into account the system properties. For instance one could exploit sparsity 
as in [WKKM06] by considering the variable sparsity correlation pattern of the 


polynomials {T l ,i 
1 ,..., njn} and k. 




€ X, j = 1, 


= 1 ,...,n 0 },{rf,j = 


5 Benchmarks 


Here, we perform some numerical experiments while solving Problem (13 ) (given 
in Section W) on several examples. In Section 5.1 we verify that the program 
of Example l] satisfies some boundedness property. We also provide examples 
involving higher dimensional cases. Then, Section [5.2| focuses on other properties, 
such as checking that the set of variable values avoids an unsafe region. Numerical 
experiments are performed on an Intel Core i5 CPU (2.40 GHz) with Yalmip 
being interfaced with the SDP solver Mosek. 


5.1 Checking boundedness of the set of variables values 

Example f. Following Example[l] we consider the constrained piecewise discrete¬ 
time dynamical system S = (X ln , X°, {A 1 , X 2 }, {T 1 , T 2 }) with X ln = [0.9,l.l]x 
[0,0.2], A' 0 = {x £ K 2 | r°(x) < 0} with r° : x hA —1, X 1 = {x £ R 2 | r 1 (x) < 0} 
with r 1 : x 1 —>■ ||x|| 2 — 1, X 2 = {x £ K 2 | r 2 (x) < 0} with r 2 = — r 1 
and T 1 : (xi,x 2 ) i-t (cnx 2 + ci 2 X 2 ,c 2 ixf + c 22 x|), T 2 : (xi,x 2 ) i-A (dnxf + 
di 2 X 2 , d 2 ix 2 + d 22 x|). We are interested in showing that the boundedness prop¬ 
erty 'P ||.||2 a holds for some positive a. 

Here we illustrate the method by instantiating the program of Example [l] with 
the following input: a\ = 0.9, a 2 = 1.1, b 1 = 0, & 2 = 0 . 2 , cn = ci 2 = c 2 i = c 22 = 











(a) m = 3 


(b) m = 4 


(c) m = 5 


Fig. 2. A hierarchy of sublevel sets P m for Example [4] 


1, dn = 0.5, di 2 = 0.4, d 2 i = —0.6 and d 2 2 = 0.3. We represent the possible 
initial values taken by the program variables (*i,x 2 ) by picking uniformly N 
points {x^\x^) {i = 1,..., N) inside the box X m = [0.9,1.1] x [0,0.2] (see the 
corresponding square of dots on Figure [2j. The other dots are obtained after 
successive updates of each point (x ( /' ) , x^) by the program of Example [I] The 
sets of dots in Figure [2] are obtained with N = 100 and six successive iterations. 

At step m = 3, Program (131 yields a solution {p 3 ,w 3 ) £ Rg[x] x R together 
with SOS certificates, which guarantee the boundedness property, that is a; £ 

=> x £ P 3 := {p 3 (x) < 0} C V\\. 111^3 => ||a :||2 < w 3 . One has p 3 {x) := 
—2.510902467—0.0050*1 —0.0148*2+3.0998*5—0.8037*5 —3.0297x?+2.5924*5 + 
1.5266*1*2 - 1.9133*5* 2 - 1.8122*1*5 + 1.6042*? + 0.0512x?x 2 - 4.4430*5*5 - 
1.8926*1*5 + 0.5464*? - 0.2084*? + 0.5866*?* 2 + 2.2410*f*i + 1.5714*5*5 - 
0.0890*1*5 ~ 0.9656*5 + 0.0098*? — 0.0320*?* 2 — 0.0232*?*| + 0.2660*?*5 + 
0.7746*5*5+0.9200*1*2+0.6411*5 (for the sake of conciseness, we do not display 
Pi and p 5 ). 

Figure [2] displays in light gray outer approximations of the set of possible 
values Xi taken by the program of Example [4] as follows: (a) the degree six 
sublevel set P 3 , (b) the degree eight sublevel set P 4 and (c) the degree ten 
sublevel set P 3 . The outer approximation P 3 is coarse as it contains the box 
[—1.5,1.5] 2 . However, solving Problem (131 at higher steps yields tighter outer 
approximations of 91 together with more precise bounds uq and w 3 (see the 
corresponding row in Table |2|. 

We also succeeded to certify that the same property holds for higher di¬ 
mensional programs, described in Example [5] (d = 3) and Example U(d = 4). 


Example 5. Here we consider X m = [0.9,1.1] x [0,0.2] 2 , r° : * 1 —> —1, r 1 : 
* 1 ^ ||*||5 — 1, r 2 = — r 1 , T 1 : (* 1 ,* 2 , * 3 ) 1 —> 1/4(0.8*5 + 1.4* 2 — 0.5x§, 1.3*i + 
0.5*5,1.4*2+0.8*|), T 2 : (*i; *2) * 3 ) | —t 1/4(0.5*i+0.4*|,—0.6*5+0.3*|, 0.5*3+ 
0.4*5) and n : * i-» ||*|||. 

Example 6. Here we consider X ln = [0.9,1.1] x [0,0.2] 3 , r° : * 1 — Y —1, 7' 1 : * 1 —»• 
||*||5 — 1, r 2 = — r 1 , T 1 : (* 1 , * 2 , *3, *4) 1 —> 0.25(0.8*5 + 1.4* 2 — 0.5*§, 1.3*i + 









0.5, x\ — O. 8 X 4 ,0.8cc§ + 1 . 4 x 4 , 1 - 3 x 3 + 0.5x1), T 2 : (xi, X 2 , X 3 , X 4 ) 1 —> 0.25(0.5xi + 
0.4x|, —0.6xf + 0.3x|, 0 . 5 x 3 + 0.4x1, —O. 6 X 3 + 0.3X 2 ) an d « : x 1 —>■ ||x|||. 


Table [l] reports several data obtained while solving Problem © at step m, 
(2 < m < 5), either for Example |4j Example [ 5 ] or Example [6] Each instance of 
Problem © is recast as a SDP program, involving a total number of “Nb. vars” 
SDP variables, with a SDP matrix of size “Mat. size”. We indicate the CPU time 
required to compute the optimal solution of each SDP program with Mosek. 

The symbol ” means that the corresponding SOS program could not be 
solved within one day of computation. These benchmarks illustrate the com¬ 
putational considerations mentioned in Section as it takes more CPU time 
to analyze higher dimensional programs. Note that it is not possible to solve 
Problem (131 at step 5 for Example [6] A possible workaround to limit this com¬ 
putational blow-up would be to exploit the sparsity of the system. 


Table 1. Comparison of timing results for Example [4] [5] and [6] 


Degree 2 m 

4 

6 

8 

10 

Example [ 4 ] 

Nb. vars 

1513 

5740 

15705 

35212 

Mat. size 

368 

802 

1404 

2174 

(d = 2) 

Time 

0.82 s 

1.35 s 

4.00 s 

9.86 s 

Example [H] 

Nb. vars 

2115 

11950 

46461 

141612 

Mat. size 

628 

1860 

4132 

7764 

CO 

II 

Time 

0.84 s 

2.98 s 

21.4s 

109 s 

Example [6] 

Nb. vars 

7202 

65306 

18480 

_ 

Mat. size 

1670 

6622 

373057 

— 

(d = 4) 

Time 

2.85 s 

57.3 s 

1534 s 

- 


Table 2. Hierarchies of bounds obtained for various properties 


Benchmark 

K 

W 2 W 3 Wa W 5 

Example 

Example 

4 

7 


x h* 0.25- ||x + 0.5||| 

639 17.4 2.44 2.02 

0.25 0.249 0.0993 -0.0777 

Example 

8 


xh-> \\T\x)-T 2 {x)\\l 

10.2 2.84 2.84 2.84 

5.66 2.81 2.78 2.78 


5.2 Other properties 

Here we consider the program given in Example [7] One is interested in showing 
that the set Xi of possible values taken by the variables of this program does 
not meet the ball B of center (—0.5, —0.5) and radius 0.5. 


















Example 7. Let consider the PPS S — (A ln , X°, {X 1 , X 2 }, {T 1 , T 2 }) with X ln = 
[0.5,0.7] x [0.5, 0.7], A' 0 = {x £ R 2 | r°(x) < 0} with r° : x —1, X 1 = {x £ 
R 2 | r 1 (x) < 0} with r 1 : x >->• ||x|| 2 — 1, A 2 = {x £ R 2 | r 2 (x) < 0} with 
r 2 = —r 1 and T 1 : (xi,X 2 ) i-t (x 2 + x|,xf + x 2 ), T 2 : (x,y) >->• (0.5xf + 
0.4x|, —0.6x 2 + 0.3x|). With k : (xi,X 2 ) >->• 0.25 — (xi + 0.5) 2 — (X 2 + 0.5) 2 , one 
has B := {x £ R 2 | 0 < k(x)}. Here, one shall prove x € 93 => k(x) < 0 while 
computing some negative a such that 93 C 'P KjCt . Note that k is not a norm, by 
contrast with the previous examples. 


At step m = 3 (resp.m = 4), Program yields a nonnegative solution wj, 
(resp. W 4 ). Hence, it does not allow to certify that 93 fl B is empty. This is 
illustrated in both Figure [3] (a) and Figure [3] (b), where the light grey region 
does not avoid the ball B. However, solving Program (131 at step m = 5 yields 
a negative bound W 5 together with a certificate that 93 avoids the ball B (see 
Figure |3] (c)). The corresponding values of w m (to = 3,4,5) are given in Table [ 2 ] 



(a) m = 3 (b) m = 4 (c) to = 5 

Fig. 3. A hierarchy of sublevel sets P m for Example [T] 


Finally, one analyzes the program given in Example [8j 


Example 8. (adapted from Example 3 in (AJ13] 1 

Let S be the PPS (AT in , AT 0 , {A 1 , AT 2 }, {T 1 , T 2 }) with A in = [-1,1] x [-1,1], 
A' 0 = {x £ K 2 | r°(x) < 0} with r° : x h> —1, A 1 = {x £ R 2 | r x (x) < 0} with 
r 1 :xh>x 2 -Xi, A' 2 = {x £ R 2 | r 2 (x) < 0} with r 2 = — r 1 and T 1 : (xi,x 2 ) '-t 
(0.687x1+0.558x2 —0.0001*xiX2, — 0.292xi + 0.773x2), T 2 : (x,y) <—> (0.369xi + 
0.532x2—O.OOOlxf, —1.27x1+0.12x2—O.OOOIX 1 X 2 ). We consider the boundedness 
property K\ := || • ||| as well as K 2 (x) := ||T 1 (x) — T 2 (x)|||. The function k 2 can 
be viewed as the absolute error made by updating the variable x after a possibly 
“wrong” branching. Such behaviors could occur while computing wrong values 
for the conditionals (e.g. r 1 ) using floating-point arithmetics. Table [ 2 ] indicates 
the hierarchy of bounds obtained after solving Problem (13 ) with to = 3,4, 5, for 


both properties. The bound w$ = 2.84 (for «i) implies that the set of reachable 
values may not be included in the initial set A ln . A valid upper bound of the 
error function H 2 is given by w$ = 2.78. 










6 Templates bases 


We finally present further use of the set P defined at Equation ([7]). This sub- 
level set can be viewed as a template abstraction, following from the definition 
in |AGGllt . with a fixed template basis p and an associated 0 bound. This rep¬ 
resentation allows to develop a policy iteration algorithm [ AGM15j to obtain 
more precise inductive invariants. 

We now give some simple method to complete this template basis to improve 


the precision of the bound w found with Problem (13). 


Proposition 2 (Template basis completions). Let ( p,w ) be a solution of 
Problem (131 and Q be a finite subset o/R[x] such that for all q £ Q, p—q £ £[x\. 
Then SR C {i g p{x) < 0, q{x) <0, \/q £ Q} C 
{a: £ | p{x) < 0, q(x) <0, Vg £ Q} is an inductive invariant. 


V C V 

' K,W — ' K,a 


and 


Proof. Let Q be the set {x £ | p(x) < 0, q(x) <0, Vq £ Q}. It is obvious 

that Q C P = {x £ M. d | p(x) < 0} and hence Q C V K}W . Now let us prove that 
Q is an inductive invariant. We have to prove that Q satisfies Equation |8| that 
is: (i) For all x £ A ln , q{x) < 0; (ii) For all i £ I, for all x £ Q fl X‘ fl A' 0 , 
q{T l {x)) < 0. For all q £ Q, we denote by ip q the element of £[x\ such that 
p — q = if q . Let us show (i) and let x £ X ln . We have q(x) = p{x) — if q (x) and 
since if q £ A[cc], we obtain, q(x) < p{x). Now from Proposition [l] and Lemmajl] 
and since (p, w) is a solution of Problem (131, we conclude that q[ x) < p(x) < 0. 


Now let us prove (ii) and let i £ 1 and x £ Q fl X 1 fl A 0 . We get q{T l {x)) = 
p(T l (x)) — ipqiT 1, (x)) and since £ i7[a;], we obtain q(T l (x)) < p(T t (x)). Using 
the fact that (p, w) is a solution of Problem (131 and using Proposition [T] and 
Lemma [l] we obtain q(T l {x)) < p{T l {x)) < p{x). Since x£QCP={y£R d \ 
p(y) < 0}, we conclude that q(T l (x)) < 0. 


Actually, we can weaken the hypothesis of Proposition [2] to construct an in¬ 
ductive invariant. Indeed, after the computation of p following Problem (13), it 
suffices to take a polynomial q such that p — q> 0. Nevertheless, we cannot com¬ 
pute easily such a polynomial q. By using the hypothesis p — q £ A[;c], we can 
compute q by sum-of-squares. Proposition [2] allows to define a simple method to 
construct a basic semialgebraic inductive invariant set. Then the polynomials de¬ 
scribing this basic semialgebraic set defines a new templates basis and this basic 
semialgebraic set can be used as initialisation of the policy iteration algorithm 
developed in [AGM15] . Note that the link between the templates generation and 
the initialisation of policy iteration has been addressed in |Adj 14|. 


Example 9. Let us consider the property 'P ||.||2 00 and let (p, w) be a solution of 
Pro blem (|13| . We have k(x) = Y^i<j<k x ^j anc ^ w + p — k = if where if £ A[.t]. 
In |RJGF12j , the templates basis used to compute bounds on the reachable 
values set consists in the square variables plus a Lyapunov function. Let us 
prove that, in our setting, Q = {x% — w, k = 1,..., d} can complete {p} in the 
sense of Proposition [ 2 ] Let k £ {1,... ,d} and let x £ R d , p(x) — (x% — w) = 
p{x) - k(x) +w+ x ) = V’(z) + Xj e £[x\. 

















7 Related works and conclusion 


Roux et al. [R.TGF12] provide an automatic method to compute floating-point 
certified Lyapunov functions of perturbed affine loop body updates. They use 
Lyapunov functions with squares of coordinate functions as quadratic invariants 
in case of single loop programs written in affine arithmetic. In the context of 
hybrid systems, certified inductive invariants can be computed by using SOS 
approximations of parametric polynomial optimization problems [LWYZ14] , In 
EH!, the authors develop a SOS-based methodology to certify that the trajec¬ 
tories of hybrid systems avoid an unsafe region. 

In the context of static analysis for semialgebraic programs, the approach 
developed in |Cou05j focuses on inferring valid loop/conditional invariants for 
semialgebraic program^ This approach relaxes an invariant generation prob¬ 
lem into the resolution of nonlinear matrix inequalities, handled with semidef- 
inite programming. Our method bears similarities with this approach but we 
generate a hierarchy of invariants (of increasing degree) with respect to target 
polynomial properties and restrict ourselves to linear matrix inequality formula¬ 
tions. In |BRCZ05j . invariants are given by polynomial inequalities (of bounded 
degree) but the method relies on a reduction to linear inequalities (the polyhedra 
domain). Template polyhedra domains allow to analyze reachability for polyno¬ 
mial systems: in jSTDG12) . the authors propose a method that computes linear 
templates to improve the accuracy of reachable set approximations, whereas the 
procedure in |DT12j relies on Bernstein polynomials and linear programming, 
with linear templates being fixed in advance. Bernstein polynomials also appear 
in [RG13] as polynomial templates but they are not generated automatically. 
In [SG09], the authors use SMT-based techniques to automatically generate 
templates which are defined as formulas built with arbitrary logical structures 
and predicate conjunctions. Other reductions to systems of polynomial equalities 
(by contrast with polynomial inequalities, as we consider here) were studied in 
|MQS04IRCK07] and more recently in |CJJK14] , 

In this paper, we give a formal framework to relate the invariant generation 
problem to the property to prove on analyzed program. We proposed a prac¬ 
tical method to compute such invariants in the case of polynomial arithmetic 
using sums-of-squares programming. This method is able to handle non trivial 
examples, as illustrated through the numerical experiments. Topics of further 
investigation include refining the invariant bounds generated for a specific sub- 
level property, by applying the policy iteration algorithm. Such a refinement 
would be of particular interest if one can not decide whether the set of variable 
values avoids an unsafe region when the bound of the corresponding sums-of- 
squares program is not accurate enough. For the case of boundedness property, 
it would allow to decrease the value of the bounds on the variables. Finally, 
our method could be generalized to a larger class of programs, involving semi¬ 
algebraic or transcendental assignments, while applying the same polynomial 
reduction techniques as in [;MAGW14| . 

c This approach also handles semialgebraic program termination 
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